import numpy as np
from scipy import optimize as op

# 随机生成10个点
m = 50
np.random.seed(10) # 固定随机数种子
randPoint = np.array([])

def randomPoint3D():
    return np.array([np.random.randint(0, 100), np.random.randint(0, 100), np.random.randint(0, 100)])

for i in range(m):
    randPoint = np.append(randPoint, randomPoint3D())
randPoint.shape = m, 3

randPoint
# 计算每个点的p

def grad_f(point):
    return np.array([2 * point[0] , 4 * point[1], 6 * point[2]]) / 2

pointP = np.array([])

for i in range(m):
    pointP = np.append(pointP, grad_f(randPoint[i]))

pointP.shape = m,3
pointP
# 计算矩阵A
def getA(m):
    A = np.zeros((m*(m-1), m))
    line = 0
    for i in range(m):
        for j in range(m):
            if i == j:
                continue
            A[line][i] = -1
            A[line][j] = 1
            line += 1
    return A

A = getA(m)
print(A)
print(A.shape)
# 计算矩阵b
b = np.array([])

for i in range(m):
    for j in range(m):
        if i == j:
            continue
        d1 = np.sqrt(np.sum((randPoint[i] - pointP[j])**2))
        d2 = np.sqrt(np.sum((randPoint[i] - pointP[i])**2))
        b = np.append(b, d1 - d2)

print(b)
b.shape
# 计算矩阵c
c = np.array([])
for i in range(m):
    c = np.append(c, 1)

print(c)
c.shape
A_ub = A
b_ub = b
res = op.linprog(c, A_ub=A_ub, b_ub=b_ub)
print(res)
